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This paper explores the use of a discrete singular convolution algorithm as a unified approach for 
numerical integration of the Fokker-Planck equation. The unified features of the discrete singular 
convolution algorithm are discussed. It is demonstrated that different implementations of the present 
algorithm, such as global, local, Galerkin, collocation, and finite difference, can be deduced from a 
single starting point. Three benchmark stochastic systems, the repulsive Wong process, the Black- 
Scholes equation and a genuine nonlinear model, are employed to illustrate the robustness and to 
test accuracy of the present approach for the solution of the Fokker-Planck equation via a time- 
dependent method. An additional example, the incompressible Euler equation, is used to further 
validate the present approach for more difficult problems. Numerical results indicate that the present 
unified approach is robust and accurate for solving the Fokker-Planck equation. 
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Much research has been done in the exploration of accurate and stable computational methods for the numerical 
solution of the Fokker-Planck equation |]]-[2^]. A detailed comparison of several different approaches was given by 
"2^ (see Ref. for a literature review). In fact, the solution of the Fokker-Planck equation, 
inear form of this equation, is still a non-trivial problem. In somewhat a broader sense, finding 
numerical solutions for partial differential equations (PDEs) is still a challenge owing to the presence of possible 
singularities and/or homoclinic manifolds that induce sharp transitions in the solutions p7[ |. These phenomena can 
be observed in many real systems such as black holes in astronomy, shock waves in compressible fluid flow, vortex 
sheets in incompressible flow associated with a high Reynolds number, and burst events in the turbulent boundary 
layer. The difficulties associated with these phenomena can often be characterized by sharp changes occurring in 
a very small spatial region which can strongly influence the global properties of the system. The presence of these 
phenomena can be extremely sensitive to numerical algorithms and can easily lead to numerically induced spatial 
and/or temporal chaos p^ . At present, there are two major classes of numerical methods for solving PDEs, namely, 
global methods and local methods. In global methods, unknown functions and their derivatives are expanded in terms 
of a finite basis set with each element having a global support. The expansion coefficients are often determined by 
the method of tau, or Galerkin, or collocation, or others. In the Tau method, the residual for a truncated expansion 
is required to be orthogonal to a subset of basis functions used in the expansion, which, together with the boundary 
conditions, determines the expansion coefficients. In the global Galerkin method, a new set of basis functions is 
constructed by the superposition of the original basis functions. The requirement of the residual be orthogonal to the 
new set of basis functions, together with the boundary conditions, determines the expansion coefficients. In the global 
collocation approach, the residual vanishes at a subset of node points of the highest order basis function used in the 
expansion. The global collocation is also called pseudospectral method. Three most important local approaches are 
finite difference, finite volume and finite element methods. In finite difference methods, the solution is interpolated 
in terms of a set of grid values; the spatial derivatives are usually approximated by algebraic expressions involving 
nearest neighbor grid points. In finite volume approaches, the emphasis is on a set of integro-differential equations 
and their associated surface and volume integrations. The values on the boundary of each "numerical molecule" are 
usually interpolated by low order schemes. The spatial derivatives are approximated in the same way as those used 
in the finite difference methods. Finite element methods form one of the most versatile classes of numerical methods. 
Depending on the system under study, finite element methods can be formulated either in terms of the method of 
weight residuals or in terms of variational principles. Usually, PDEs are integrated by using a set of trial functions, 
each with a small region of support. The solution is represented by linear superpositions of these trial functions. 

Global methods are highly localized in their spectral space, but are unlocalized in the coordinate space. By contrast, 
local methods have high spatial localization, but are delocalized in their spectral space. In general, global methods are 
much more accurate than local methods, while the major advantage of local methods is their flexibility for handling 
complex geometries and boundary conditions. Moreover, the use of global methods is usually restricted to structured 
grids, whereas, local methods can be implemented to block-structured grids and even unstructured grids. 

There were hectic debates among the numerical computation communities over the advantages and disadvantages 
of various numerical methods in the past a few decades. These debates stimulated the development of powerful 
numerical methods for a wide variety of science and engineering applications. Such development has, in association 
with the availability of inexpensive high-performance computers, led to the establishment of numerical simulations as 
an alternative approach for researches and applications. The connection of various numerical methods has always been 
an important research topic. Finlayson discussed the relation between the Galerkin and the Ritz variational principle 
p9| . Canuto et al rearranged their spectral basis functions so that some global collocation method can be regarded as a 
special case of certain global Galerkin methods ||3(|[. F ornberg addressed the common feature between pseudospectral 
methods and high order finite difference methods [|31|. The connection between global and local methods can also be 
realized in the framework of the method of weighted residual by choosing trial functions of cither piecewise Lagrange 
polynomials or global Lagrange polynomials. The connection of methods of finite element, finite difference and finite 
volume is now well understood [ p2| . However, to our knowledge, none has reported a unified scheme for the discussion 
of all of the abovementioned methods. 

In previous work | |25[ |, we proposed a discrete singular convolution (DSC) algorithm and demonstrated its use for 
the numerical solution of Fokker-Planck equation via eigenfunction expansions. The DSC algorithm was shown to 
be a potential numerical approach for Hilbert transform, Abel transform. Radon transform and delta transform. 
Three standard problems, the Lorentz Fokker-Planck equation, the bistable model and the Henon-Heiles system, were 
utilized to test the accuracy, reliability, and speed of convergence of the DSC-eigenfunction approach. All results 
were in excellent agreement with those of previous methods in the field. Recently, the DSC algorithm has been 
successfully tested for integrating the sine-Gordon equation with initial values close to homoclinic orbits [^3|, which is 
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extremely difficult to compute because of the possible presence of numerical chaos . Excellent results are obtained 
for solving the Navier-Stokes equation and for engineering structural analysis . The purpose of the present paper 
is twofold. First, we study the unified features of the DSC algorithm for treating partial differential equations. This 
is accomplished by focusing on the DSC kernels of the delta type and their approximations. Second, we explore the 
use of the DSC as a unified approach for solving the Fokker-Planck equation via direct explicit time propagations. 
The eigenfunction expansion approach provides a Schrodinger-equation type picture for the understanding of the 
Fokker-Planck equation. However, its use is restricted to a certain class of Fokker-Planck operators (essentially for 
the Fokker-Planck operators their equivalent Schrodinger potentials are bounded from below). The present direct 
approach is applicable to a wider class of problems. These two DSC-based approaches have the same level of accuracy 
for the numerical solution of the Fokker-Planck equation. They are complementary to each other for solving a wide 
variety of Fokker-Planck systems arising from practical situations. 

This paper is organized as the follow. The unified features of the DSC algorithm arc discussed in Section II. We 
demonstrate that, the present DSC algorithm provides a unified framework for solving the Fokker-Planck equation, 
and partial differential equations in general. In particular, we show that various different implementations of the DSC 
algorithm, such as global, local, Galerkin, collocation, and finite difference, can be deduced from a single starting 
point. The application of the present DSC approach to the solution of the Fokker-Planck equation and Euler equation 
is presented in Section III. We use four examples to illustrate the present approach. The first example is the repulsive 
Wong process which is useful for testing the ability of handling monomodality-bimodality transition. The second 
example is the Black-Scholes equation for option derivatives. This is interesting stochastic model for option pricing in 
financial market. The third case treated is a nonlinear stochastic model which has certain connection to a mean-field 
model for self-organization processes in biological systems such as muscle contraction. Notably, all of these problems 
are treated by an explicit time-propagation approach in contrast to the eigenfunction expansion used in our previous 
work p^ ]. Since the abovementioned examples are of strong parabolic type, we consider an additional problem, the 
incompressible Euler equation, to further validate the DSC approach for more difficult problems. The incompressible 
Euler equation is chosen because its equations for velocity vector and pressure field are of strong hyperbolic type and 
elliptic type, respectively. Thus, this last example is complimentary to other three examples from the point of view 
of numerical analysis. This paper ends with a discussion. 

II. PROPERTIES OF THE DISCRETE SINGULAR CONVOLUTION 

This section presents the properties of the discrete singular convolution (DSC) algorithm for solving differential 
equations. The first subsection addresses the unified features of the DSC algorithm in the line of the method of 
weighted residuals. Relevant properties of DSC trial functions are discussed the second subsection. 

A. Unified features 

Without the loss of generality, it is assumed that at a fixed time, a stochastic process is governed by a differential 
equation. To solve the differential equation, one can start with either by approximating the original differential 
operator or by approximating the actual solution of the differential equation while maintaining the original differential 
operator. The latter is accomplished by explicitly defining a functional form for approximations. Let us assume that 
the differential equation has the form 

Cu{x) = fix), x£n, (1) 

where >C is a linear operator and u{x) is the unknown solution of interest. Here f{x) is a known force term, Q denotes 
the domain over which the differential equation applies. 

The approximate solution is sought from a finite set of N DSC trial functions of a given resolution a, denoted by 
with M being the half width of support of each element. Here cr is a regularization parameter for improving the 
regularity of the set. The case of regularization free is easily obtained by setting a ^ oo. Elements of the set S^'^'^ 
can be explicitly given by {<l>^a-i' 'f'a'a-N}- For a given computational domain, the resolution parameter a 

is determined by N. 

An important property of the DSC trial functions {4>aa-k} that when the trial function is free of regularization, 
each member of the set is a reproducing kernel at highest resolution 

lim <<,^fe,77>=7?(xfe), (2) 
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where < •, • > denotes the standard inner product. In fact, if an appropriate basis is used for (f> and the hmit on a is 
taken, of each resolution can be a reproducing kernel for functions bandlimited to appropriate sense. In general, 
we require the low pass filter property that for given a ^ O.a ^ and M ^ 

<^i'„.,k,V>~vM. (3) 

This converges uniformly when the resolution is refined, e.g., a — > oo. Many examples of such DSC trial functions are 
given in Refs. and |^6|. Further discussion on these functions is given in the next subsection. Equations (||) and 
are special requirements satisfied by the DSC kernels of delta type . 
In the present DSC approach, an approximation to the function of interest u{x) can be expressed as a linear 
combination 

N 

UZa'\^)^Y.^..o;k^'la;k{^), (4) 

where x is an independent variable and Ua,a-k is a DSC approximation to the solution wanted at point Xk- This 
structure is due to the DSC trial function property and it dramatically simplifies the solution procedure in 
practical computations. 

In this formulation, we choose the set S^'^^ a priori, and then determine the coefhcients {Ua,a;k} so that U^'^'^{x) 
is a good approximation to u{x). To determine Ua.a-k, we minimize the amount by which U^'^^{x) fails to satisfy 
the original governing equation (|l]). A measure of this failure can be defined as 

Rff{x)^CUZf{x)~f{xl (5) 

where R^'^^{x) is the residual for particular choices of resolution, regularization and half width of the support. Note 
that Eq. is similar to the usual statement in the method of weighted residuals. However, the approximation 
U^]^^{x) is constructed by using the DSC trial functions, fp^a-kix), in the present treatment. Let Eq. (|l|) and its 
associated boundary conditions be well-posed, then there exists a unique solution u{x) which generally resides in an 
infinite-dimensional space. Since the DSC approximation U^'^^ is constructed from a finite-dimensional set, it is 
generally the case that U^'^'^{x) ^ u{x) and therefore R^'^'^{x) ^ 0. 

Galerkin. We seek to optimize _R^'*^(x) by forcing it to zero in a weighted average sense over the domain fl. A 
convenient starting point is the Galerkin 

RZfix)<p^':,^,4x)dx^o, K':,.,4x) e s^:f , (6) 

n 

where the weight set S^, can be simply chosen being identical to the DSC trial function set S^'^ . We refer Eq. 
as a DSC-Galerkin statement. 
Collocation. First, we note that in view of Eq. (|^), the present DSC-Galerkin statement reduces to a collocation 
one at the limit of a' 



hm / <f {x)K'^^,.^i{x)dx = <f {xi) = 0, (7) 

where {xi} is the set of collocation points. However, in digital computations, we cannot take the above limits. It 
follows from the low pass filter property of the DSC trial functions, Eq. (0), that 



RZf{x)K'''..'A^)d^ « RZfi^i) « 0- (8) 



It can be proven that for appropriate choice of S^, ''^f , the first approximation of Eq. (H) converges uniformly. The 
difference between the true DSC-coUocation, 

lim Rl-^\xi)^Q, (9) 



and the Galerkin induced collocation, diminishes to zero for appropriate DSC trial functions. 

Global and local. Global approximations to a function and its derivatives are realized typically by a set of truncated 
i^(a,6) function expansions. It is called global because the values of a function and its derivatives at a particular 
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point Xi in the coordinate space involve the full set of grid points in a computational domain fl. Whereas a local 
method does so by requiring only a few nearest neighbor points. In the present DSC approach, since the choices of M 
and/or M' are independent of N, one can choose M and/or M' so that a function and its derivatives at a particular 
point xi are approximated either by the full set of grid points in the computational domain or just by a few nearest 
neighbor grid points. In fact, this freedom for the selection of M endows the DSC algorithm with controllable accuracy 
for solving differential equations and the flexibility for handling complex geometries. 

Finite Difference. In the finite difference method, the differential operator is approximated by difference operations. 
In the present approach, the DSC-coUocation expression of Eq. (||) is equivalent to a 2M + 1 (or 2M) term finite 
difference method. This follows from the fact that the DSC approximation to the nth order derivative of a function 
can be rewritten as 



'dxfl 



k+M 

•.=k-M 



u{xi), 



(10) 



where c^^ are a set of DSC weights for the finite difference approximation and are given by 

(11) 

Obviously, for each different choice of (p^^^., we have a different DSC-finite difference approximation. Hence, the 
present DSC approach is a generalized finite difference method. This DSC-finite difference was tested in previous 
studies When M — 1, the DSC-finite difference approximation reaches its low order limit and the resulting 

matrix is tridiagonal. In this case, the present DSC weights c^^ can always be made exactly the same as those 

of the second order central difference scheme (i.e. 2S'' '^j ~ ^^^^ order derivative and -jr , — , -jr for the 

second order derivative. Here A is the grid spacing.) of the standard finite difference method by appropriately 
choosing the parameter a. However, even in this case, the DSC-finite difference approximation does not have to be 
the same as the standard finite difference scheme and can be optimized in a practical application by varying a. 



B. DSC trial functions 



There are many DSC trial functions that satisfy Eq. (^. The requirement of Eq. (^ can be regarded as an 
approximate reproducing kernel or quasi reproducing kernel. The reason for using an approximate reproducing kernel 
can be understood from the following analysis of Shannon's kernel '^'"^"^^ _ Shannon's kernel is a delta sequence 

lim^^=5(x), (12) 

a^oo TTX 

where S{x) is the delta distribution which can be regarded as a universal reproducing kernel because its Fourier 
transform is the unit. However, such a universal reproducing kernel cannot be directly used in digital computations 
because it is a distribution (Precisely, it is belongs to Sobolev space of order -1, .) and it does not have a value 
anywhere in the coordinate space. Therefore, in certain sense, constructing a reproducing kernel in an appropriate 
L^{a,b) space is equivalent to finding a sequence of approximation of the delta distribution in the L'^{a,b). In fact. 
Shannon's kernel is an element of the Paley- Wiener reproducing kernel Hilbert space 

m - r fivf^^^^^dy, VfeBl, (13) 

j—oo ""i-^ y) 

where V/ G B"^ indicates that, in its Fourier representation, the function / vanishes outside the interval [— tt, tt]. 
What is important for digital computations is the fact that the Paley- Wiener reproducing kernel Hilbert space has a 
sampling basis Sk{x) 

sin7r(a; - yk) i, w, ^ ^ (,a\ 

Sk[x) = — p yk = fc, vfc e Z, (14) 

7r(a; - yk) 

where symbol Z denotes the set of all integers. Expression (|l^) provides a discrete representation of every (continuous) 
function in B^ 
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fix) = E f(yk)Sk{x), V/ G Bl (15) 

This is Shannon's samphng theorem and is particularly important to information theory and the theory of sampling. 
Note that Shannon's kernel is obviously interpolative on Z 

(16) 

where Sn.m is the Kronecker delta function. Computationally, being interpolative is of particular importance for 
numerical accuracy and simplicity. 

In wavelet analysis, '^'"^^^^ is Shannon's wavelet scaling function and its Fourier transform is a characteristic 
function, i.e. it is an unsmoothed, ideal low pass filter. In physical language, it is a projection to the frequency 
subband [— tt, tt]. By the Heisenberg uncertainty principle, such a (sharp) projection must be an infinite impulse 
response (IIR) filter. The usefulness of such a filter is limited because it is de- localized in the coordinate space 
and requires infinitely many sampling data. In practical computations, a truncation is required, which leads to large 
truncation error and even worse, numerical instability. To improve the smoothness and regularity of Shannon's kernel, 
we introduce a regularization 

$.(x) = ^-^^^R^ix) (a > 0), (17) 



TTX 



where Ra- is a regularizer which has properties 



and 



lim R^{x) = 1 (18) 



RaiO) = 1. (19) 



Here Eq. (|8|) is a general condition that a regularizer must satisfy, while Eq. (|T^) is specifically for a delta regularizer, 
which is used in regularizing a delta kernel. Various delta regularizers can be used for numerical computations. An 
excellent one is the Gaussian 



Rcr{x) ~ exp 



2(72 



(20) 



An immediate benefit of the regularized Shannon's kernel, Eq. (p^, is that its Fourier transform is infinitely differ- 
entiable because the Gaussian is an element of the Schwartz class functions. Qualitatively, all kernels of the Dirichlet 
type oscillate in the coordinate representation. Specifically, Shannon's kernel has a long tail which is proportional 
to i, whereas, the regularized kernels decay exponentially fast, especially when the a is very small. In the Fourier 
representation, regularized Shannon's kernels have an "optimal" shape in their frequency responses. Of course, they 
all reduce to Shannon's low pass filter at the limit 

^ , , sin TTX sin TTX ,„ , 

lim $o-(a;) = hm e ^ = . (21) 

(T^oo (T^oo TTX TTX 



Quantitatively, one can examine the normalization of $cr(a;) 

J <i>„{x)dx = 4,(0) 



oo / T\k / \2fc 



erf 



^^^k\(2k + l) Vx/2 



71 



= 1 - \ e — e ^ ""^dt 

\ TT a Jo 

^l-erfc(^) (23) 
1, (24) 
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where erf(z) — e * dt is the error function and erfc(z) is the complementary error function. Note that for 

a given a > 0, erfc(^) is positive definite. Thus, $ct(0) is always less than unity except at the limit of ct — > oo. 
Therefore, ^a{x) is no longer a reproducing kernel. However, we argue that ^aix) is an approximate reproducing 
kernel because when we choose a ^ v^/tt, which is the case in many practical apphcations, the residue term, erfc(^), 

approaches zero very quickly. As a result, $(t(0) is extremely close to unity. As trial functions, regularized Shannon's 
kernels do not form a sampling basis. They are no longer orthogonal in general. However, they just slightly miss the 
orthogonality and the requirement of a basis. 

For numerical computations, it turns out that the approximate reproducing kernel has much less truncation errors 
for interpolation and numerical differentiations. Qian and the present author | |3^ have recently given the following 
theorem for truncation errors. 

Theorem Let / be a function / e L'^(R) H C*(i?) and bandlimited to B, {B < ^, A is the grid spacing). For a fixed 
t E R and ct > 0, denote g{x) — f{x)Hk{^^^), where Hk{x) is the fcth order Hermitc polynomial. If g{x) satisfies 



g'{x) < g{x) 

for x>t+{Mi- 1)A, and 

g'{x) > g{x) 

hr x <t - M2A, where Mi,M2 G TV, then for any s G Z+ 



{x-t) 



t) 



(25) 



(26) 



< a/3 



rii+Mi 



sin^(i-nA) ^ (t-nAf, 
fit-nA) 



L27ra(f -S)exp(^^^^^) 



l/WII 



" ^ Via 



exp 



(MiA)2 



) 



ll/Wlh 



exp(i^) 



-M2A 

via- 



(27) 



where superscript, (s), denotes the sth order derivative. The proof and detailed discussion (including a comparison 
with the truncation errors of Shannon's sampling theorem) are given in Ref. and are beyond the scope of this 
paper. 

This theorem provides a guide to the choice of M, a and A. For example, in the case of interpolation (s = 0), if 
the L2 norm error is set to 10^'' (77 > 0), the following relations can be deduced from Eq. (p7|) 



and 



r(7r - BA) > a/4.6177, 



M , 

— > V4-61?7, 
r 



(28) 



(29) 



where r — a j A (The choice of a is always proportional to A so that the width of the Gaussian envelope varies with 
the central frequency). The first inequality states that for a given grid size A, a large r is required for approximating 
high frequency component of an function. The second inequality indicates that if one chooses the ratio r = 3, then 
the half bandwidth M ^ 30 can be used to ensure the highest accuracy in a double precision computation (77 — 15). 
However, for lower accuracy requirement, a much smaller half bandwidth can be used. In general, the value of r is 
proportional to M . The use of M values is determined by the accuracy requirement. This theoretical estimation is 
in excellent agreement with a previous numerical test ||3q|. 
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III. ILLUSTRATIVE CALCULATIONS 



In this section, we illustrate the use of the present approach for solving the Fokker-Planck equation and the 
incompressible Euler equation. Many DSC kernels are discussed in the previous work |25|j3^ ] and they can be used as 
the DSC trial functions. For simplicity, we focus on three DSC kernels, a regularized Shannon's kernel (RSK), 



.M ( ^ sinf(a;-a;fc) 



a regularized Dirichlet kernel (RDK), 
and a regularized Lagrange kernel (RLK) 



sin [-J (a; - Xk) 



{x - Xkf 
2(72 



(2m + 1) sin 



A 2m+l 



exp 



(x - Xkf 



2cr2 



(30) 



(31) 



2m 



Xk 



■ exp 



(a: - XkY 



2cr2 



(32) 



for our numerical test. Note that the resolution is given by a = -J which is the frequency bound in the Fourier 
representation. The goal of this section is to test the present method for the solutions of the Fokker-Planck equation 
via time propagation and the incompressible Euler equation. For the numerical solution of the Fokker-Planck equation, 
we choose a = 3.8A for the RSK and RDK, a = 2.8A for the RLK, with tt/A being the resolution. In fact, a wide 
range of a values can be used to deliver excellent results. The half bandwidth, M, can be chosen to interplay between 
the local limit and the global limit and is set to 40 in all calculations. Finally, m controls the order of the regularized 
Dirichlet and Lagrange kernels and is set to 40 in all calculations (note that the selection of m is independent of the 
grid used in the computation). It is noted that all of the abovementioned DSC trial functions are of Schwartz class 
and are capable of auto-regularizing when used as integral kernels. The fourth order explicit Runge-Kutta scheme is 
used for time discretization. Details of these computations are described in the first three subsections. For treating 
the incompressible Euler equation, many other DSC parameters are tested as indicated in the last subsection. Double 
precision is used in all calculations. 



A. The repulsive Wong process 

One of important stochastic systems is the repulsive Wong process p8|- ^ , p0| , given by 

dx = 27tanh(a;)dt + V2dFt, 
where dFt is the Gaussian white noise which has the standard statistical properties 

< dFt >= 

and 

< dFt,dFr >^Si\t-T\). 



(33) 
(34) 
(35) 



The repulsive Wong process is Markovian due to the deriving Gaussian white noise term. Its transition probability 
density is governed by the Fokker-Planck equation of the form Pq-ptI 



df{x,t) „ d[t&nh{x) f{x,t)] , d^f{x,t) 
- -27- 



with the usual initial condition 



dt dx 



/(x,0) = S{x - xo), 



dx'^ 



(36) 



(37) 



and the normalization 
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f{x,t)dt = l. (38) 



For 7 = 1, the solution [ p9y40| of the Fokker-Planck equation ( p6D is analyticaUy given by (for xq — 0) 



e 'I* + e 



(39) 



where x± = ±2t are centers of two moving Gaussians. Here the superposition of two Gaussians gives rise to a 
monomodahty-bimodahty transition as time increases. The Wong process is useful for illustrating the connection 
between stochastic processes and quantum measurements. It is also useful for distinguishing spectrum differences 
between the Master equation and its Fokker-Planck equation approximations. 

The accurate simulation of the Wong process is not a simple task because of the monomodality-bimodality transition. 
Two Gaussian peaks centered at x± = ±2t move apart as time increases. The computational domain is to be 
sufhciently large in order to avoid boundary reflection (Otherwise, more complicated techniques, such as absorption 
boundaries, are to be implemented.). In the present computations, the resolution is chosen as -^^ = IOtt. The initial 
functions are approximated by a unit impulse function located at 0. The equation ( ^6| ) is integrated up to 10 time 
units with a time increment of 0.001. The errors for a wide range of propagation times are listed in TABLE I and are 
measured by error norms of L2 and Loo • It is seen that the present unified approach is extremely accurate and reaches 
machine precision. All of the DSC kernels behave very similar to each other and provide the same level of accuracy 
and speed of convergence. In fact, other DSC kernels, such as regularized modified Dirichlet kernel, provide similar 
results. The results of the RSK and RDK are slightly more accurate than those of the RLK. It is evident that the 
present unified method, in associated with the DSC trial functions, is capable of delivering extremely high accuracy 
and numerical stability for the Wong process. To our knowledge, the DSC solution for this system is the best to the 
date. 

B. The Black-Scholes equation 

The Fokker-Planck equation and stochastic analysis have interesting applications in mathematical modeling of 
financial market option pricing. Consider a writer of a European call option on a stock, he is exposed to the risk of 
unlimited liability if the stock price rises acutely above the strike price. To protect his short position in the option, 
he should consider purchasing certain amount of stock so that the loss in the short position in the option is offset by 
the long position in the stock. In this way, he is adopting a hedging procedure. A hedge position combines an option 
with its underlying asset so as to achieve the goal that either the stock protects the option against loss or the option 
protects the stock against loss. This risk-monitoring strategy has been commonly used by practitioners in financial 
markets. The most well-known stochastic model for the equilibrium condition between the expected return on the 
option, the expected return on the stock and the riskless interest rate is the Black-Scholes equation [|l| 

where S is the asset price which undergoes geometric Brownian motion, c{S, t) the call price, v the volatility and r 
the constant riskless interest rate. Black-Scholes equation is a fundamental equation in finance and economics and is 
also an excellent example application of stochastic analysis. By a simple transformation 

x^l-aS, (41) 

and 

f{x,t) = e^'c{x,t), (42) 
the Black-Scholes equation is transformed into the Fokker-Planck equation of the standard form 

The numerical simulation of the Black-Scholes equation and its generalized versions is an important issue in financial 
analysis and computational finance community p^-p5[|. Essentially, all existing numerical methods are tested for 
potential usefulness in estimating the option derivatives because both computational accuracy and efficiency are very 
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important to option modeling and risk estimation. In the present time-dependent approach, the resolution is set to 

2 

■g- = 27r and the time increment is chosen as 0.01. For simplicity, we choose ^ = 0.5 and r = 0.7 in our calculations. 
We chosen our initial distribution as a unit impulse function located at a; = 0, which is a poor approximation to the 
true delta distribution. Obviously, had one started with a smooth initial function, or used a denser grid, one would 
have obtained much higher accuracy at earlier times as well. We have verified this computationally, but these results 
are not presented. Both L2 and Loo error analyses are used to evaluate the quality of the DSC approach, the results 
of which are listed in TABLE II. To our knowledge, the present time-dependent DSC approach provides the most 
accurate numerical results yet obtained for the Black-Scholes equation. 

As in the first example, three DSC kernels provide extremely similar results in solving the Black-Scholes equation. 
This is not an isolated coincidence. In fact, we can come up a number of other DSC kernels with all of their results 
being very similar to those of the present three kernels. 



C. A nonlinear stochastic model 



To illustrate the accuracy and robustness of the present approach further, we choose the following nonlinear stochas- 
tic model 

df{x,t) ^ d[{Lox + 6<x{t)>)f{x,t)] d^f{x,t) 

dt dx dx^ ' ^ ' 

where < x{t) > is the first moment of the distribution function 



< x{t) >^ / xf{x,t)dx, (45) 

J —00 

and UJ, 6 and D are constant. The initial probability distribution is also given by 

f{x,0)^5{x-xo). (46) 

Equation ( ^ ) is a true nonlinear stochastic model since the instantaneous position average depends on the distribution 
function. This is one of few analytically soluble nonlinear systems which are very valuable for testing new numerical 
approaches. For example, Drozdov and Morillo have recently employed this system to test their K-point Stirling 
interpolation formula finite difference method p3| . The exact solution to Eq. (p[) is 



fi^' = J— 77^ exp 



{x- < x{t) >Y 
2v{t) 



(47) 



where < x{t) > and v[t) are analytically given by 

< x{t) >= xoe-('^+®)* (48) 

and 

v{i) = - (1 - e-2-*) , (49) 

UJ 

respectively. Obviously, v{t) is the theoretical value of the second moment M2{t) 

M2{t)^<x^{t)> - <x{t)>^, (50) 

which can also be used as a measure of computational accuracy. 

In the present computations, the resolution is chosen as ^ — ^tt. The time increment is taken as At = 0.005. In 
this example, the errors are measured by error norms of Li and Loo from which all other error norms, such as the 
L2 error norm, can be interpolated. The Li and Loo errors are listed in TABLE III, for D = 0.1, = 1,9 = 2 and 
xq = 2.0422. The initial accuracy of computations is hindered by the poor approximation of the impulse function 
to the Dirac delta function. However, the auto-regularization property of the Schwartz class trial functions enables 
the numerical integration to stabilize at smooth solution and eventually reach the machine precision at a slightly late 
time. 
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D. The Euler equation 



All cases considered in the last three subsections are of strong parabolic type with a solution which becomes more 
and more flat and smooth as time increases. In this subsection, we consider an additional problem, the incompressible 
Euler equation, to confirm that the results obtained for the Fokker-Planck equation are not due to the parabolic nature. 
We also use this example to demonstrate the inter-connection between the collocation and the finite difference, and 
between the local and the global. It is hoped that this additional example helps to build confidence for using the 
DSC approach for treating more difficult problems. Conceptually and numerically, it is convenient to describe the 
incompressible Euler equation from the point of view of the incompressible Navier-Stokes equation 

Vv = Vp+— V^v, (51) 
ot Re 

V-v==0, (52) 

where v is the velocity field vector, p is the pressure field and Re is the Reynolds number. The Euler equation is 
attained by setting Re= oo. Finding a general solution to the Euler equation is not an easy job. In the present study, 
we consider a solution domain of [0, 2tt] x [0, 2tt] with periodic boundary conditions. Under such a constraint, the 
Navier-Stokes equation (|5|) exists an exact solution 

u{x, y,t) — ~ cos(x) sin(y)e~ 
v{x, y, t) = sin(.T) cos(y)e~ 

1 _ 
p{x,y,t) ^ ~-[cos{2x) +cos{2y)]e , (53) 

where (u, v) are the velocity components in the x-direction and ?/-direction, respectively. Note that, for the Euler 
equation, the solution ( |53| ) does not decay with time. 

In the present study, we use a standard approach for treating the incompressible Navier-Stokes equation, i.e. 
deriving a Poisson equation for the pressure from the incompressible condition. The velocity fields are iterated by 
using the implicit Euler scheme. At time tn+i, there are two coupled equations for the velocity fields 



(54) 



i^V^-i^).-^=prV2 + ^«, (55) 



and a Poisson equation for the pressure 
Here 5", Sy and Sp are given by 
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S:^-— + iuV^+v-u;) 

5;' = -£ + (u"<+«X) 

s; = ^« + <) - - K)' - (57) 

At each time i„+i, the pressure field p»+i/2 is 

solved according Eq. ( p6| ) from the known velocity field vector (m", tj"). 
The velocity field vector (w"+^, 1;"+^) is then updated according to Eqs. (^J) and (|5^). These linear algebraic equations 
are solved by using a standard (LU decomposition) solver. 

The derivatives in Eqs. (^-56) are computed by using the generalized finite difference scheme, Eq. (10), and the 



required finite difference weights are given by Eq. (^_1|). The involved trial functions, (b'^ „^ are given by the regularized 
Shannon's kernel (RSK) [Eq. (|30|)], the regularized Dirichlet kernel (RDK) [Eq. (|3l|)], and the regularized Lagrange 
kernel (RLK) [Eq. (^]. Here m = 40 is used for both RDK and RLK. We choose a small time increment (At = 0.001) 
so that the main error is caused by the spatial discretization. The number of grid points in each dimension is chosen 
as TV = 4, 8, 16 and 32 in various test calculations. The a value is specified as a = 4- = — For a given iV, 

the matrix half bandwidth, Af, can be chosen as M < N. In particular, if M = N, the approach has a global (full) 
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computational matrix. For all M < N , the matrix is banded. In the present DSC approach, the connection between 
the global and the local can be easily achieved by selecting an M value for a given TV. In particular, \i M <^ N , the 
DSC approach behaves truly like a finite difference scheme. To achieve optimal (or near optimal) accuracy, the a is 
chosen in proportional to M and A. When M — 32, 16, 8, 4, 2 and 1, ^ are chosen as 3.2, 2.5, 1.8, 1.2, 0.9 and 0.6 for 
both RSK and RDK, and 2.8, 2.0, 1.6 1.0 0.8 and 0.6 for RLK. We compute the L2 and Loo errors of u for a number 
of combinations of N and M and the results are listed in TABLE IV for 4 different times {t = 0.5, 1.0, 1.5,2.0). A 
good consistent in accuracy among solutions at different times (or equivalently, over 2000 iterations) is observed. The 
DSC results are quite accurate when N = M = A and are of machine precision when N = M = 32. It is interesting to 
note that for fixed M = 4, the results of = 32 and A^ = 4 differ little in accuracy. We also checked the DSC-finite 
difference approximation at the tridiagonal matrix limit {Ad = 1) and the result is very good for A^ = 4 (i.e. a total 
of 4 interior points in a (27r)^ box). 

IV. DISCUSSION 

The main purpose of this paper is to discuss the unified features of the discrete singular convolution (DSC) algorithm 
p5| . It is found that the implementations of the DSC algorithm into a number of computational methods can be 
deduced from a single starting point, the method of weighted residual. This chain of deduction provides a imified 
approach for solving the Fokker-Planck equation and other differential equations in general. Some of these deduction 
relations are novel to our knowledge. 

We demonstrate that by adjusting the support of the DSC trial functions, the DSC algorithm can be easily im- 
plemented either as a local method or as a global method. For this reason, the present DSC approach has global 
method's accuracy while maintains local method's flexibility for handling complex boundary and geometry. In fact, 
the solution of the Fokker-Planck equation of a previous paper j2^ was obtained by using the global limit. Whereas, 
in the present computations, a local approximation is used for all Fokker-Planck problems. A comparison between 
global and local DSC treatments is given in solving the Euler equation. 

We also show that the DSC implementations of Galerkin and collocation are computationally equivalent, i. e. the 
collocation, Eq. (||) can be deduced from the Galerkin, Eq. (||) because of the choice of the DSC trial functions. 
Galerkin methods have a profound influence to the theory of approximations. Both spectral methods and finite 
element methods are often formulated in the framework of the Galerkin approach. There has been a great deal of 
argument about advantage and disadvantage of the Galerkin in comparison to many other methods. The present 
DSC approach might provide a unified framework for the discussion of these methods. 

The present Galerkin-induced collocation scheme provides a nature base for the realization of finite difference 
methods. High order finite difference is not a new idea in numerical approximations |3|]. However, the mathematical 
constructions of high order finite difference schemes often become too cumbersome to use in practical applications as 
the order increases. The present DSC approach provides a simple, systematic algorithm for the generation of finite 
difference schemes of an arbitrary order. The implementation of this finite difference is demonstrated in solving the 
Euler equation with a number of different matrix bandwidths. 

Recently, wavelet theory and techniques have had great success in signal processing, data compression, and telecom- 
munication. Two most important features of the wavelet theory are multiresolution analysis and time- frequency lo- 
calization. Their potential applications in solving partial differential equations have been extensively explored |46H5l|] 
in hope to come up with unifled approaches for numerical approximations. However, before wavelet approaches can 
be of practical use, a number of technical difficulties are to be overcome. In our view, the first difficulty is the imple- 
mentation of boundary conditions in a multiresolution setting. The second difficulty is the requirement of sufficiently 
high wavelet regularity to provide sufficiently weak solutions. Moreover, there is a lack of general and systematic 
numerical algorithms for incorporating wavelets in an efficient manner. Nevertheless, the wavelet multiresolution 
analysis still has great potential for developing adaptive grid and multigrid algorithms. The present DSC algorithm 
is closely related to the wavelet theory [^,^. In fact, the DSC kernels have a feature in common with wavelets in 
terms of time-frequency (position-momentum) localization. However, unlike in a wavelet algorithm, multiresolution 
analysis is feasible but it is not required in the DSC algorithm. 

In contrast to our earlier work dealing with the application of the DSC approach to the Fokker-Planck equation via 
an eigenfunction expansion |25| , we have explored in this paper a DSC-based time-propagation approach for solving 
the Fokker-Planck equation. Three typical DSC kernels, the regularized Shannon's kernel (RSK), the regularized 
Dirichlet kernel (RDK), and the regularized Lagrange kernel (RLK), are used as trial functions in the framework of 
the present method. Four benchmark examples are chosen to demonstrate the usefulness and to test the accuracy 
of the present DSC approach. The first example is the repulsive Wong process. This is used for objectively testing 
the ability of handling the monomodality-bimodality transition. The Wong process requires a large computational 
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domain to ensure that the boundary reflection of density flux has httle influence in a highly accurate computation. 
By using reasonable resolution, regularization and a quite large time increment, the present approach performs very 
well in characterizing the transition. In fact, the present unified approach delivers machine accuracy at an early time. 

The Black-Scholes equation of option pricing was chosen as the second numerical example. This financial equation 
can be regarded as a reaction-diffusion equation, although, its derivation was entirely based stochastic analysis. By 
using a simple; transformation, the Black-Scholes equation is converted into the standard form of the Fokker-Planck 
equation which admits an analytical solution. The present numerical results for the Black-Scholes equation are 
obtained by using three different DSC kernels with a reasonable resolution and relatively large time mesh. The 
extremely high accuracy in the present calculation indicates that the DSC-based unified algorithm is a valuable 
potential approach for various option pricing simulations. 

The third example treated is a nonlinear stochastic model. The effective potential of the corresponding Fokker- 
Planck equation is time depended through the first order moment of the transition probability density. Despite of 
the nonlinearity and poor approximation of the initial density distribution, the numerical solutions quickly settle to 
a smooth, stable and correct distribution after a few iterations. This is due to the fact that the DSC trial functions 
are chosen as Schwartz class functions and they are capable of auto-regularizing when used as integration kernels. 
Our results are of machine precision at a relatively late time. To our knowledge, this is the best numerical solution 
to this nonlinear Fokker-Planck equation to the date. These illustrative calculations indicate that the present unified 
approach is extremely accurate, efficient and robust for numerical simulations of stochastic systems. 

A common feature in the abovementioned Fokker-Planck equation is that the equation is of strong parabolic type 
and the solution decays as time increases. Therefore, it is necessary to employ an additional example to validate the 
present DSC algorithm further for handling more complicated partial differential equations. To this end, we choose the 
incompressible Euler equation with its velocity field equations being of strong hyperbolic type and a derived equation 
for the pressure being of elliptic type. A standard implicit Euler scheme is used for the time discretization and at 
each time tn+i, linear algebraic equations are constructed by using the collocation method. In the present approach, 
carrying out differentiations in the collocation is equivalent to implementing the finite difference weights computed 
from the DSC trial functions. We test the DSC algorithm by using 4, 8, 16 and 32 grid points in each dimension 
in association with many different half matrix bandwidths (M= 4, 8, 16 and 32). As expected, the DSC algorithm 
achieves its highest accuracy at the global limit (M = TV) for each given N. The machine precision is reached when 
N = M = 32. Very good results are also obtained for many banded matrix calculations. We believe that the feature 
of being able to provide both global and local approximations in one formulation is of practical importance for large 
scale computations. 

Although this paper emphasizes the connection of a few computational methods and the unified features of the 
DSC approach, it claims neither that all computational methods are the same, nor that the DSC algorithm engulfs 
all methods. For example, it is still not clear whether the DSC algorithm is applicable in adaptive and unstructured 
grids (progress is made on a DSC-multigrid method). The reader is urged to keep the distinction of various methods 
in mind and maintain a perspective. 
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TABLE I. Errors for solving the repulsive Wong process 



Time 


RSK 


RDK 


RLK 


L2 




L2 




L2 




0.10 


1.94(-09) 


2.32(-09) 


1.94(-09) 


2.32(-09) 


1.94(-09) 


2.32(-09) 


0.25 


4.67(-ll) 


3.73(-ll) 


4.67(-ll) 


3.73(-ll) 


4.67(-ll) 


3.72(-ll) 


0.50 


3.47(-12) 


2.05(-12) 


3.47(-12) 


2.05(-12) 


3.47(-12) 


2.03(-12) 


0.75 


8.76(-13) 


4.89(-13) 


8.76(-13) 


4.95(-13) 


8.83(-13) 


4.60(-13) 


1.00 


3.57(-13) 


1.83(-13) 


3.57(-13) 


1.91(-13) 


3.78(-13) 


1.91(-13) 


2.00 


6.29(-14) 


2.69(-14) 


5.92(-14) 


3.00(-14) 


2.18(-13) 


8.92(-14) 


3.00 


4.72(-14) 


2.04(-14) 


3.60(-14) 


1.63(-14) 


2.86(-13) 


9.76(-14) 


4.00 


5.23(-14) 


1.97(-14) 


3.76(-14) 


1.39(-14) 


3.53(-13) 


1.12(-13) 


5.00 


4.44(-14) 


1.19(-14) 


3.16(-14) 


1.12(-14) 


4.10(-13) 


1.21(-13) 


6.00 


5.87(-14) 


1.88(-14) 


5.23(-14) 


1.64(-14) 


4.68(-13) 


1.32(-13) 


7.00 


7.63(-14) 


2.38(-14) 


7.30(-14) 


2.16(-14) 


5.25(-13) 


1.43(-13) 


8.00 


9.27(-14) 


2.78(-14) 


9.08(-14) 


2.58(-14) 


5.80(-13) 


1.54(-13) 


9.00 


7.19(-14) 


1.99(-14) 


6.24(-14) 


1.79(-14) 


6.33(-13) 


1.63(-13) 


10.0 


6.98(-14) 


1.70(-14) 


5.25(-14) 


1.39(-14) 


6.87(-13) 


1.73(-13) 
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TABLE II. Errors for the numerical solution of the Black-Scholes equation 



Time 


RSK 


RDK 


RLK 


L2 




L2 




L2 




1 


1.85(-03) 


1.14(-03) 


1.83(-03) 


1.12(-03) 


2.18(-03) 


1.41(-03) 


2 


1.19(-04) 


8.55(-05) 


1.18(-04) 


8.45(-05) 


1.42(-04) 


1.06(-04) 


3 


4.57(-06) 


2.89(-06) 


4.48(-06) 


2.84(-06) 


6.44(-06) 


4.29(-06) 


4 


2.53(-07) 


1.68(-07) 


2.46(-07) 


1.63(-07) 


4.23(-07) 


2.84(-07) 


5 


1.81(-08) 


1.15(-08) 


1.75(-08) 


1.10(-08) 


3.66(-08) 


2.43(-08) 


6 


1.63(-09) 


1.09(-09) 


1.56(-09) 


1.04(-09) 


4.02(-09) 


2.74(-09) 


7 


1.81(-10) 


1.26(-10) 


1.71(-10) 


1.19(-10) 


5.48(-10) 


3.92(-10) 


8 


2.41(-11) 


1.59(-11) 


2.25(-ll) 


1.48(-11) 


9.03(-ll) 


6.05(-ll) 


9 


3.83(-12) 


2.55(-12) 


3.54(-12) 


2.37(-12) 


1.76(-11) 


1.15(-11) 


10 


7.82(-13) 


6.25(-13) 


7.30(-13) 


5.85(-13) 


4.01(-12) 


2.84(-12) 


20 


4.86(-14) 


2.70(-14) 


4.91(-14) 


2.78(-14) 


4.84(-14) 


2.74(-14) 
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TABLE III. Errors for solving the nonlinear model 



Time 


RSK 


RDK 


RLK 


Li 




Li 




Li 




1 


4.13(-01) 


3.09(-02) 


6.46(-01) 


4.71(-02) 


4.67(-02) 


3.98(-03) 


2 


4.66(-01) 


3.68(-02) 


7.41(-01) 


5.86(-02) 


2.48(-03) 


2.27(-04) 


3 


2.03(-01) 


1.68(-02) 


3.26(-01) 


2.68(-02) 


7.24(-04) 


6.51(-05) 


4 


4.88(-02) 


4.92(-03) 


7.97(-02) 


8.04(-03) 


1.23(-04) 


1.31(-05) 


5 


7.41(-03) 


7.84(-04) 


1.22(-02) 


1.29(-03) 


1.68(-05) 


1.82(-06) 


6 


1.02(-03) 


1.12(-04) 


1.68(-03) 


1.84(-04) 


2.26(-06) 


2.48(-07) 


7 


1.38(-04) 


1.52(-05) 


2.28(-04) 


2.50(-05) 


3.06(-07) 


3.36(-08) 


8 


1.87(-05) 


2.05(-06) 


3.09(-05) 


3.38(-06) 


4.14(-08) 


4.54(-09) 


9 


2.53(-06) 


2.77(-07) 


4.17(-06) 


4.58(-07) 


5.59(-09) 


6.13(-10) 


10 


3.42(-07) 


3.75(-08) 


5.65(-07) 


6.19(-08) 


7.56(-10) 


8.29(-ll) 


20 


3.64(-14) 


3.11(-15) 


4.75(-14) 


4.88(-15) 


1.00(-13) 


1.51(-14) 
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TABLE IV. Errors for solving the Euler equation 









RSK 


RDK 


RLK 


iV 


M 


Time 


Li 




Li 




Li 




4 


1 


0.5 


3.16(-2) 


6.12(-2) 


3.15(-2) 


6.12(-2) 


3.09(-2) 


6.01(-2) 






1.0 


3.10(-2) 


6.02(-3) 


3.09(-2) 


6.01(-3) 


3.03(-2) 


5.88(-3) 






1.5 


3.06(-2) 


5.90(-3) 


3.05(-2) 


5.89(-3) 


2.99(-2) 


5.73(-3) 






2.0 


3.04(-2) 


5.77(-3) 


3.04(-2) 


5.76(-3) 


2.99(-2) 


5.57(-3) 




2 


0.5 


1.27(-2) 


2.48(-3) 


1.27(-2) 


2.48(-3) 


1.29(-2) 


2.44(-3) 






1.0 


3.10(-2) 


6.02(-3) 


1.31(-2) 


2.66(-3) 


1.29(-2) 


2.54(-3) 






1.5 


1.37(-2) 


2.86(-3) 


1.37(-2) 


2.87(-3) 


1.32(-2) 


2.66(-3) 






2.0 


1.44(-2) 


3.08(-3) 


1.45(-2) 


3.08(-3) 


1.35(-2) 


2.78(-3) 




4 


0.5 


9.33(-3) 


1.70(-3) 


9.32(-3) 


1.69(-3) 


9.34(-3) 


1.70(-3) 






1.0 


9.43(-3) 


1.79(-3) 


9.42(-3) 


1.79(-3) 


9.44(-3) 


1.79(-3) 






1.5 


9.62(-3) 


1.88(-3) 


9.61(-3) 


1.88(-3) 


9.64(-3) 


1.88(-3) 






2.0 


9.92(-3) 


1.96(-3) 


9.91(-3) 


1.95(-3) 


9.93(-3) 


1.96(-3) 


8 


8 


0.5 


1.30(-4) 


4.26(-5) 


1.33(-4) 


4.36(-5) 


1.24(-4) 


3.40(-5) 






1.0 


1.52(-4) 


5.13(-5) 


1.54(-4) 


5.25(-5) 


1.54(-4) 


4.79(-5) 






1.5 


1.82(-4) 


6.10(-5) 


1.83(-4) 


6.26(-5) 


1.92(-4) 


5.66(-5) 






2.0 


2.17(-4) 


7.13(-5) 


2.16(-4) 


7.32(-5) 


2.34(-4) 


6.72(-5) 


16 


16 


0.5 


6.30(-10) 


2.37(-10) 


6.75(-10) 


2.76(-10) 


1.23(-8) 


3.63(-9) 






1.0 


6.76(-10) 


2.40(-10) 


6.82(-10) 


2.87(-10) 


1.56(-8) 


5.16(-9) 






1.5 


8.00(-10) 


2.65(-10) 


7.57(-10) 


3.18(-10) 


1.99(-8) 


6.76(-9) 






2.0 


9.68(-10) 


3.35(-10) 


8.81(-10) 


3.49(-10) 


2.48(-8) 


8.53(-9) 


32 


4 


0.5 


5.25(-4) 


2.10(-4) 


5.24(-4) 


2.10(-4) 


2.37(-3) 


7.14(-4) 






1.0 


7.40(-4) 


2.88(-4) 


7.41 (-4) 


2.90(-4) 


2.96(-3) 


1.00(-3) 






1.5 


1.04(-3) 


4.25(-4) 


1.04(-3) 


4 27f-4) 


3.75(-3) 


1.33(-3) 






2.0 


1.40(-3) 


5.87(-4) 


1.40(-3) 


5.89(-4) 


4.62(-3) 


1.67(-3) 




8 


0.5 


1.78(-6) 


7.32(-7) 


1.93(-6) 


7.95(-7) 


9.50(-7) 


3.05(-7) 






1.0 


2.41(-6) 


1.06(-6) 


2.62(-6) 


1.15(-6) 


1.24(-6) 


5.22(-7) 






1.5 


3.23(-6) 


1.39(-6) 


3.51(-6) 


1.51(-6) 


1.64(-6) 


7.73(-7) 






2.0 


4.17(-6) 


1.79(-6) 


4.52(-6) 


1.92(-6) 


2.09(-6) 


1.03(-6) 




16 


0.5 


6.95(-ll) 


2.93(-ll) 


9.52(-ll) 


4.00(-ll) 


1.90(-10) 


5.73(-ll) 






1.0 


7.48(-ll) 


3.23(-ll) 


1.03(-10) 


4.42(-ll) 


2.40(-10) 


8.15(-11) 






1.5 


8.06(-ll) 


3.51(-11) 


l.ll(-lO) 


4.81(-11) 


3.04(-10) 


1.08(-10) 






2.0 


8.67(-ll) 


3.79(-ll) 


1.19(-10) 


5.19(-11) 


3.76(-10) 


1.36(-10) 




32 


0.5 


1.02 (-14) 


6.99(-15) 


1.36(-14) 


1.21(-14) 


1.10(-14) 


7.88(-15) 






1.0 


2.03(-14) 


1.45(-14) 


2.51(-14) 


1.80(-14) 


2.22(-14) 


1.54(-14) 






1.5 


2.98(-14) 


1.88(-14) 


3.68(-14) 


2.74(-14) 


3.54(-14) 


2.22(-14) 






2.0 


4.05(-14) 


2.31(-14) 


5.04(-14) 


3.13(-14) 


4.80(-14) 


3.12(-14) 
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